Probing ergodicity in granular matter 
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When a granular system is tapped, its volume changes. Here, using a well-defined macroscopic 
protocol, we prepare an ensemble of granular systems and track the statistics of volume changes as a 
function of the number of taps. This is in contrast to previous studies, which have focused on single 
trajectories and assumed ergodicity. We devise a new method to assess the convergence properties 
of a sequence of ensemble volume histograms and introduce a reasonable approximate version of 
an invariant histogram. We then compare these invariant histograms with histograms generated by 
sampling a long trajectory for one system and observe nonergodicity, which we quantify. Finally, 
we use the overlapping histogram method to assess potential compatibility with Edwards' canonical 
assumption. Our histograms are incompatible with this assumption. 
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I. INTRODUCTION 

Controlling the packing density of a vibrated granu- 
lar system is of industrial and scientific interest. In this 
paper, we focus on sequences of packings that are in me- 
chanical equilibrium. For such packings, the effect of 
shaking depends solely on the tapping stage of the shak- 
ing protocol. This stage is uniquely parameterized by 
the extra momentum given to the system that sets the 
amount of dilation allowed after each tap [l] . 

The protocol used to prepare the granular system is 
also of primary importance. Currently, no macroscopic 
preparation protocol exists that allows one to choose the 
initial microscopic configuration exactly. Due to the in- 
evitable dispersion in the initial conditions, the only re- 
producible quantity that one can measure in experiments 
is the probability of observing a given volume after a 
given number of taps. 

During the last twenty years, properties of the time 
averages of packing volumes and their fluctuations have 
been studied extensively [T-6 . Recently McNamara et 
al. [7 carried out a careful analysis of volume histograms 
sampled over time. However, if the systems under study 
are not ergodic — an a priori assumption in most of these 
studies — then time averages over fluctuations for one sys- 
tem need not be the same as averages over the initial con- 
ditions. To tackle this problem, one should ideally study 
an ensemble of systems prepared in the same macroscopic 
way and follow the evolution of their volume with time. 
To our knowledge, such a study is still lacking. 

In this paper, we use numerical simulations to study 
the evolution of an ensemble of granular systems pre- 
pared with the same macroscopic protocol. We analyse 
the rate of convergence of these histograms and find that, 
under certain conditions, individual trajectories are non- 
ergodic. Finally, we find that our results are incompatible 
with Edwards' hypothesis, namely, that the dependence 
of the volume statistics on preparation protocol can be 
captured by a simple Boltzmann weight [8j. 



II. "EQUILIBRIUM" DISTRIBUTION 



A. Numerical simulation 



The shaking experiment is simulated using a dissi- 
pative, event-driven MD simulation of frictionless hard 
spheres subject to gravity. We use lateral periodic bound- 
ary conditions and place a wall at the bottom of the sim- 
ulation box. 

A known problem with event-driven schemes applied to 
dissipative systems is that a linear restitution coefficient 
makes the simulation "freeze" locally in higher-density 
regions [9]. To avoid this problem, we use a nonlinear 
restitution coefficient, so that the spheres do not loose 
all their kinetic energy, even after many collisions. When 
the system relaxes after an external perturbation (e.g. a 
single tap), the final state will have a well-defined struc- 
ture, but the spheres will still move locally — albeit with 
a typical displacement amplitude that is much smaller 
than their size. We assume that the packings thus ob- 
tained are representative of the typical stable configura- 
tions that are obtained in real experiments at sufficiently 
high densities. 

We start every simulation by generating an equili- 
brated fluid configuration of 1000 hard spheres at pack- 
ing fraction ^ = 0.2 in a cubic simulation box. We then 
switch on gravity and let the system reach mechanical 
equilibrium, as explained above. This is the preparation 
stage of the simulation. Once the packing has settled, we 
"tap" the system by accelerating the bottom layer of the 
spheres impulsively. Concretely, we add to these spheres 
the vertical velocity As/ga^ where A is a dimensionless 
tapping amplitude, g is the acceleration of gravity, and a 
is the diameter of the spheres. We repeat this process at 
regular time intervals, which are longer than the typical 
time it takes the packing to settle after a perturbation. 
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B. Volume histogram evolution 

To study the ensemble statistics of our system, we pre- 
pare about 50000 stable initial conditions as described 
above. The corresponding distribution of the initial vol- 
umes is denoted po{V). For a given amplitude A, the 
volume histogram after i tapping steps is denoted (V). 
The volume of a packing is defined as the smallest axis- 
aligned cuboidal volume that completely contains all the 
spheres, and whose bottom face lies in the plane of the 
wall. 

In analogy with the Gibbs approach in statistical me- 
chanics, we define the statistical "equilibrium" ensemble 
distribution as the asymptotic volume distribution ob- 
tained when the number of taps tends to infinity. If there 
is a well-defined time-scale for approach to the steady 
state, then this procedure should yield a good approxi- 
mation of the invariant distribution. 
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FIG. 1. Typical volume histogram evolution for an ensemble 
of about 50000 replicas of a system of 1000 spheres. The 
picture shown corresponds to an amplitude A, equal to 6. 

In Figjljwe show the evolution of a typical volume his- 
togram as a function of the number of tapping steps. 
Starting from a very broad volume distribution, the 
ensemble histograms become narrower, and eventually 
reach an invariant shape. 



C. Convergence analysis 

To quantify the convergence oberved for the volume 
histograms in Fig^ we want to assess to what ex- 
tent a sequence {p^ (^)}i=i,2,.. becomes independent of 
i for large i. For this purpose, we introduce the 2- 
sample Kolmogorov-Smirnov (KS) statistic D[pi^p2] be- 
tween two one-dimensional histograms pi and p2 [10 : 



D[pl,p2] 



sup{|Fi(F) 

V 



■F2{V)\}, 



(1) 



where Fi and F2 are the respective cumulative distri- 
bution functions associated with histograms pi and p2. 
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FIG. 2. Convergence analysis. Plot of the KS distances 

A A 

D[pi , pi^k] running from 1 to 25 being represented by the 
different colors) for different typical amplitudes correspond- 
ing to different relaxation regimes: A = 6.0 (top), A = 7.5 
(middle) and A = 8.0 (bottom). The dashed line in each plot 
corresponds to the critical value e* for histograms built from 
50000 data points. 



This "distance" is commonly used as a test for the null 
hypothesis that the two histograms pi and p2 are differ- 
ent realisations of the same underlying distribution. This 
hypothesis can be rejected with 99% of certainty if pT] : 



D[pi,p2] > 1.63 



ni + n2 
nin2 



: e*(ni,n2), 



(2) 



where ni and n2 are the respective number of samples 
used to build pi and p2- From Eq.([2|, we expect that 
the KS statistic will decrease as the accuracy of each 
histogram increases (i.e., e*(ni, 77.2) decreases), unless the 
two histograms sample different underlying distributions. 

In practice, we seek the existence of an equilibrium 
step, denoted by /(A), past which the deviation from 
the invariant histogram is within the statistical noise. In 
Fig 2I we plot the KS distances D[p^ , p^j^^ as a function 
of both i (x-axis) and k (color code) . We also display the 
critical value e* = 0.0109 corresponding to histograms 
built from 50000 samples (dashed line, Eq.([2|). We iden- 
tify I {A) as the first step for which all colored points 
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(from black for /c = 1 to bright yellow for k = 25) are be- 
low the dashed line. We emphasize that it is not enough 

A A 

for D[p- , p^_^i] to be less than e*. For example, as shown 
in Fig|2j in the case of A = 7.5, this weaker condition is 
satisfied (black points) as early as the 30th tap, whereas 
our convergence criterion is fulfilled only after the step 
1(7.5) = 164. 

We also notice that the relaxation time increases when 
going from A = 6.0 to A = 7.5, whereas it decreases 
when going from A = 7.5 to A = 8.0. 



relaxation time. These two limits can only be reconciled 
if the relaxation time depends nonmonotonically on tap- 
ping amplitude, as shown in Figs. [2]and[3j 



III. TESTS ON THE INVARIANT 
HISTOGRAMS 

A. Ergodicity analysis 
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FIG. 3. Invariant properties. Red curve (triangles): relax- 
ation number of steps I{A) vs tapping amplitude, A. The 
line is a cubic spline interpolation, and is meant to guide the 
eye. Blue curve (circles): invariant average volume, (V), vs 
A. 

In Figjsj we plot the mean invariant volumes, (V), and 
their corresponding equilibrium step, /(A), for the am- 
plitudes we have tested. We notice that (V) increases 
with tapping amplitude, as found in previous studies 
[11 [31 |5H7]. The values of I (A) are consistent with the 
number of steps required for the ensemble average vol- 
ume to reach a plateau value for a given amplitude A 
(data not shown). As stressed previously, we note that it 
is not an always decreasing function of A. To our knowl- 
edge, this is the first time that a nonmonotonic depen- 
dence of relaxation time scales on tapping amplitude has 
been observed. Normally, one would expect a decrease 
of the relaxation time as the tapping amplitude is in- 
creased m |3l ISHZl , although a nonmonotonic dependence 
on A was already suggested in Ref.[6] but in a different 
context. Our current understanding of this dependence 
is the following. At sufficiently low A, the ensemble be- 
haves like a quenched glass. Within this regime, higher 
amplitudes result in the system exploring a larger local 
basin, which increases the time to sample the basin and, 
hence, the relaxation time. On the other hand, if we 
were to tap on the system strongly enough so that each 
shaking step is akin to our preparation protocol, then we 
should not observe any change in the ensemble histogram 
between successive taps, which would yield a vanishing 



FIG. 4. Ergodicity analysis. Plot of the distance D (n) as a 
function of the number of uncorrelated volume values, n, for 
different amplitudes, A. The plain black line corresponds to 
the threshold value e*(n) calculated from Eq.([2|. The sam- 
pling is done after having tapped 200 times on the system. 

To explain the behaviour observed in Figs. [2]and[3j we 
have proposed in the preceding section that some condi- 
tions lead to nonergodicity. In this section, we test this 
proposal quantitatively. 

Denote by p^iV) a volume histogram built from n un- 
correlated volume values taken from a very long trajec- 
tory of one system tapped with an amplitude A. In prac- 
tice, we consider a single trajectory for a system tapped 
about 10^ times and we compute the corresponding vol- 
ume correlation function. We then use the latter to select 
n uncorrelated values belonging to this trajectory. Un- 
der the ergodic hypothesis, this "time" histogram and 
the invariant ensemble histogram, Pi>/(^)C^), both sam- 
ple the same underlying distribution. To test this idea, 
we again use the KS statistic. Let D (n) be the quan- 
tity D[p^{V),p.^j^^^{V)]. If D (n) is bigger than the 
rejection value e*(n) calculated from Eq.Q, with m = n 
and 712 = 50000, then the ergodic hypothesis is rejected. 
Recall that the KS statistic decreases with increasing n 
if and only if the null hypothesis, ergodicity, is true. Oth- 
erwise, D (n) saturates at a finite value when n is large 
enough. Figkl shows that for high values of n, D (n) 
indeed tends to saturate. In that case, the difference be- 
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tween the saturation value of D (n) and the curve e* (n) 
at a fixed n serves as a useful measure of nonergodic- 
ity that we will call the ergodicity gap. Fig j4] shows that 
for tapping amplitudes A = 4.5, 6.5 and 7.5, trajectories 
are nonergodic, with ergodicty gaps that increase with 
decreasing tapping amplitude. In contrast, for A = 7.9 
and 8.0, trajectories are ergodic according to our crite- 
rion. These findings are consistent with our current un- 
derstanding of the nonmonotonic trend for I{A) observed 
in Figs. [2]and|3l 



B. Compatibility with Edwards' prior 

In 1989, Edwards et al. [8 suggested that if a granular 
material is submitted to a protocol that yields different 
packing volumes then the probability of occurence of a 
given stable configuration is proportional to a Boltzmann 
weight e~^^, where V is the actual volume of the packing 
configuration. Since then, many studies have tested this 
idea and its consequences, but no consensus has been 
reached [3 H HMZ] • 

In their study, McNamara et al. tested Edwards' prior 
by generating time-sampled volume histograms for dif- 
ferent tapping amplitudes and looking directly at ratios 
between between pairs of histograms. This procedure 
is known as the overlapping histogram method [TF: if 
the logarithm of the ratio between two histograms is 
linear, then there is a potential compatibility with Ed- 
wards' canonical assumption. Otherwise, this assump- 
tion should be rejected, at least for the tested protocol. 
McNamara et al. found a linear behaviour for the log 
ratios, for the amplitudes that they tested both numeri- 
cally and experimentally. Here, we apply the same test 
on our own invariant ensemble histograms. 

Fig|5] shows the log ratio of two invariant ensemble 
histograms {A = 7.5 with A = 8.0 for the main plot 
and A = 4.5 with A = 6.5 in the inset). We find that 
over a limited range of volumes where the histograms 
overlap significantly, there is indeed a linear decrease in 
the log ratio of the two invariant histograms, consistent 
with the findings in Ref.[7 . However, we also find that 
for high volumes, this ratio tends to saturate, whereas 
for low volumes, it rapidly decreases. These trends are 
present for all pairs of invariant volume histograms that 
we have tested, including those where the two histograms 
are apparently "closeby". 

It is worth noting that the authors of Ref. [7] have al- 
ready suggested the possibility that the observed linearity 
is only true locally. Moreover a close look at their own 
figures reveals deviations from linearity for very high and 
very low volumes that are consistent with our own obser- 
vations (Figj5|. 

Overall, Edwards' canonical hypothesis as a global 
property is not compatible with the protocol we are test- 
ing. Although this is not the first time that a strong 
disagreement with Edwards' theory has been found [19 , 
to our knowledge, this is the first time that a full his- 
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FIG. 5. Overlapping Histograms. The red points correspond 
to the plot of the logarithm of the ratios pj^j^ ^ / Pi^i^ q (main 
figure) and pj^i^ ^/pi>iQ 5 (inset). We superimposed the cor- 
responding invariant histograms in both the main figure and 
the inset. The histogram in green always corresponds to the 
smallest amplitude in a given pair. The plain lines are the 
best linear fits through the midpoints within the crossover 
regions. 



togram analysis reports incompatibility with Edwards' 
canonical assumption for vibrated granular matter. Our 
study does not rule out the possibility of a local com- 
patibility with Edwards' assumption but this is already 
different from Edwards' orignal theory. 



IV. CONCLUSION 



In this letter, we introduced an ensemble volume satis- 
tic for a simulated vertically vibrated granular system. 
To quantitatively assess the properties of the generated 
sequences of ensemble histograms, we used the KS test. 
This allowed us to devise a convergence criterion for a 
sequence of histograms. Subsequently, we tested the er- 
godicity of a tapped system as a function of the tapping 
amplitude, A, and found clear evidence for nonergodicity 
when the tapping amplitude is low. Finally, we tested the 
compatibility of our invariant histograms with Edwards' 
hypothesis and concluded that it is not compatible with 
our simulation protocol. 

We should point out that the results found in this pa- 
per depend a priori on the chosen tapping protocol and 
also on the preparation stage. Although the dependence 
of our findings on the preparation protocol is hard to 
predict, the tools that we introduced to characterize our 
ensemble histograms can be used to test any numerical 
or experimental protocol. 
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